function [bmus,errors] = neuralGasBmus(M, D, selectedBmus)

dataVectors = size(D,1);
units = size(M,1);

% selected bmus
if ((nargin < 3) || (isempty(selectedBmus))), 
  selectedBmus = 1; 
end

bmus = zeros(dataVectors,length(selectedBmus));
errors = bmus;

% set the value of blen lower. 
mask = ones(size(M,2),1);
blen = min(size(M,1),dataVectors);

W1 = (mask*ones(1,dataVectors));

% constant matrices
WD = 2*diag(mask)*D';   % constant matrix
dconst = ((D.^2)*mask); % constant term in the distances

i = 0; 
while i+1<=dataVectors, 
  % calculate distances 
  inds = [(i+1):min(dataVectors,i+blen)]; i = i+blen;      
  
  %euclidean distance
  Dist = (M.^2)*W1(:,inds) - M*WD(:,inds); % plus dconst for each sample
  
  % find the bmus and the corresponding quantization errors
  if all(selectedBmus==1), [error bmu] = min(Dist); else [error bmu] = sort(Dist); end
  if units==1, bmus(inds,:) = 1; else bmus(inds,:) = bmu(selectedBmus,:)'; end
  errors(inds,:) = error(selectedBmus,:)' + dconst(inds,ones(length(selectedBmus),1));
end  

errors = sqrt(errors);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
